// 版权所有2009 Go作者。版权所有。
// 此源代码的使用受BSD样式
// 许可证的约束，该许可证可以在许可证文件中找到。

package math

// go代码是从
// http:
// 
// ======================================================
// /太阳微系统公司1993年版权所有。保留所有权利。
// 
// 由Sun Microsystems，Inc.旗下的SunSoft公司开发。
// 使用、复制、修改和分发本
// 软件的权限是免费授予的，前提是保留本通知
// 的权利。
// ==============================================

// /Cbrt返回x的立方根。特殊情况为：
// Cbrt（±0）=±0 
// Cbrt（±Inf）=±Inf 
// Cbrt（NaN）=NaN 
func Cbrt(x float64) float64 {
	if haveArchCbrt {
		return archCbrt(x)
	}
	return cbrt(x)
}

func cbrt(x float64) float64 {
	const (
		B2             = 696219795                   // /（664-0.03306235651）*2**20 
		C              = 5.42857142857142815906e-01  // 19/35=0x3FE15F15F15F1 
		D              = -7.05306122448979611050e-01 // /-864/1225=0xBFE691DE2532C834 
		E              = 1.41428571428571436819e+00  // 99/70=0x3FF6EA0EA0EA0EA0EA0EA0EA0EA0EA0F 
		F              = 1.60714285714285720630e+00  // 45/28=0x3FF9B6DB6DB6DB6DB6E 
		G              = 3.57142857142857150787e-01  // 5/14=0x3FD6DB6DB6DB6DB6DB6DB6DB7
	)
	switch {
	case x == 0 || IsNaN(x) || IsInf(x, 0):
		return x
	}

	sign := false
	if x < 0 {
		x = -x
		sign = true
	}

	// 粗cbrt到5位
	t := Float64frombits(Float64bits(x)/3 + B1<<32)
	if x < SmallestNormal {
		// 次正常数
		t = float64(1 << 54) // 设置t=2**54 
		t *= x
		t = Float64frombits(Float64bits(t)/3 + B2<<32)
	}

	// 新cbrt到23位
	r := t * t / x
	s := C + r*t
	t *= G + F/(s+E+D/s)

	// 切到22位，使其大于cbrt（x）
	t = Float64frombits(Float64bits(t)&(0xFFFFFFFFC<<28) + 1<<30)

	// 一步牛顿迭代到53位，误差小于0.667ulps 
	s = t * t // t*t为精确
	r = x / s
	w := t + t
	r = (r - t) / (w + r) // r-s为精确
	t = t + t*r

	// 恢复符号位
	if sign {
		t = -t
	}
	return t
}
